Shaped high frequency vibratory source

ABSTRACT

The present invention is a method of processing seismic data in which one or more seismic vibrators are activated with one or more pilot signals and vibrator motions are recorded along with seismic data. Vibrator signatures are computed from measured vibrator motions, such as the ground force signal. A desired impulse response is specified from either a measured vibrator motion or from test data or field data from a location near the location from which the seismic data was acquired. A deconvolution filter is computed from the impulse response and the vibrator signature. Alternatively, a single separation and deconvolution filter is derived from the impulse response and from vibrator signatures from multiple vibrators and sweeps. The deconvolution or deconvolution and separation filter is used to process the seismic data. The vibrators are then moved to a new location, and the activation is repeated.

FIELD OF THE INVENTION

This invention relates to the field of seismic data acquisition andprocessing of vibratory source data. Specifically, this invention is amethod of optimizing vibratory source data to accurately represent theoutput that would be derived from impulse sources.

BACKGROUND OF THE INVENTION

Seismic vibrators are commonly used to produce the source signalsnecessary in the geophysical exploration for hydrocarbons. In field use,seismic vibrators are excited with a pilot signal that is typically awave train that varies in frequency, referred to as a sweep, and lastsseveral seconds. The excitation of the vibrator is typically adjusted bya feedback loop controlled by a ground force signal that is computedfrom signals measured by accelerometers on the base plate and on thereaction mass of the source.

In this application of seismic vibrators, seismograms are generated bycross correlating the data recorded at various receiver locations withthe pilot signal. This cross-correlation step compresses the impulseresponse of the data from the several seconds associated with the sweepto a few tens of milliseconds and thereby better approximates the signalthat would be recorded by an ideal impulsive source. This step isfollowed by standard seismic data processing steps, such as surfaceconsistent statistical deconvolution, static corrections, noisefiltering, bandpass filtering, and imaging.

Five categories of problems with conventional seismic vibratortechnology have long been recognized in industry. First, thecross-correlation process results in a pulse with undesirablecharacteristics, including a widened main lobe and sidelobes that appearas oscillations on either side of the main lobe. Second, the output ismixed phase, as a result of combining the correlation process, whichresults in a zero phase wavelet, with the earth attenuation filter andrecording instruments, which are minimum phase.

This leads to several problems. For example, correlated data, unlikeimpulsive-source data, do not have well-defined arrival times. Inaddition, such data are no longer causal. Other processing techniques,such as statistical or spiking deconvolution, assume that the data arecausal and minimum phase, and for that reason the processing results maynot be accurate. Third, the pilot signal used in the correlation isgenerally substantially different from the actual signal put into theground. The actual signal contains harmonics and other nonlinearitiesarising from the mechanics of the vibrator and its coupling with theground. Processing the data with the pilot signal does not allow thoseharmonics and nonlinearities to be removed, which therefore appear asnoise in the processed data. Fourth, acquiring seismic data is expensiveand a major cost is associated with the number of source stations thatcan be used. Traditional vibrator technologies record only one sourcestation at a time. Methods that allow acquisition using multiple sourcestations simultaneously would speed acquisition and reduce costs. Fifth,to increase the energy put into the ground two or more vibrators aretypically used at each source station. However, spacing limitationsresult in the multiple vibrators forming an array that can limit thehigh frequencies in the data and thereby reduce resolution. Elevationchanges may also limit the ability to correct for time shifts, alsoreferred to as static corrections, between the vibrators. As furtherdescribed in the next several paragraphs, industry has focused asubstantial amount of effort in attempting to overcome theselimitations.

The U.S. patents issued to Trantham, U.S. Pat. No. 5,400,299 andAndersen, U.S. Pat. No. 5,347,494 disclose methods that principallyaddress the first category of problems. The methods result in improvedimpulse wavelet shapes. Trantham's disclosure also provides a causal andminimum phase impulse after removal of the vibrator pilot signal.However, processing is done with the pilot, which is only anapproximation of the signature imparted into the ground. Also,Trantham's approach preferably requires pre-whitening the signal, inwhich white noise is added, to stabilize a spectral division from whichthe vibrator signature deconvolution filter is generated. This stepprevents numerical errors in the division, but also may cause phasedistortions and adds a precursor to the processed seismic data.

Andersen's method involves the choice of a sweep power spectrum thatleads to an impulse response with a simple shape and an optimum lengthafter correlation. Unlike conventional sweeps, which start around 8 Hz,Andersen's sweep starts at frequencies near 1 Hz. The presence of theselower frequencies results in a more desirable wavelet. This solution,however, only addresses the first of the problems noted above. Inaddition, the use of high-resolution wavelets requires that the sweeprate changes rapidly with time and is not realizable with standardhydraulic vibrators.

U.S. Pat. No. 5,550,786 issued to Allen discloses a method that uses ameasured accelerometer signal or signals, such as the ground force, fromeach vibrator and sweep in an inversion process instead of a correlationwith a pilot signal. The measured signals are related to the actualsignal imparted into the ground by a transfer function of minimum phase,which is obtained by the process of minimum phase statisticaldeconvolution, which is commonly used in processing land data. Stepsinclude inversion (also referred to as spectral division), bandpassfiltering, and spiking deconvolution. A model trace is processed to makea phase correction of the deconvolved data.

Allen addresses the first three shortcomings discussed above. However,Allen's method applies an inverse filter by the process of spectraldivision, and it is recognized that a problem with such inverse filtersis that at frequencies where the measured signals are small, the filterwill apply a large gain and amplify any recorded noise. If the signal iszero, then inversion will be unstable. The data can be pre-whitened byadding a small amount of constant noise to stabilize the inversion, butthe added noise can distort the phase of the data. Because processingtechniques such as spiking deconvolution assumes that input data areminimum phase, the output results may not be predictable when suchdistortion is present. Allen attempts to solve these problems by usingbandpass filters to reduce the noise outside the vibrator sweep band,and by processing a model trace in order to be able to calculate a phasecorrection. As will be understood to those skilled in the art, apreferable method would avoid the use of bandpass filters and phasecorrections by eliminating the need to pre-whiten the data.

U.S. Pat. Nos. 5,719,821 and 5,721,710 issued to Sallas et al. discuss amatrix inversion scheme to separate the outputs from individualvibrators. The number of sweeps can be equal to or greater than thenumber of vibrators. This method solves a set of linear equations, whichincludes the measured motions from each vibrator and each sweep, todetermine an optimal filter for inversion and separation. AlthoughSallas addresses the shortcomings listed above, the problems withinversion and phase errors discussed above in relation to Allen remainunsolved. U.S. Pat. No. 6,161,076 issued to Barr et al (2000) is similarto the prior work of Allen and Trantham. Barr misstates that Allen isusing a single accelerometer signal instead of the ground force signal,and claims the use of a filter to convert the data to short-durationwavelets, as did Trantham. It is understood in the art that this processis equivalent to inversion followed by a bandpass filter, and thus theproblems noted above are unaddressed. Barr specifically discloses usingthe harmonics or non-linear distortion to construct a wavelet of equalor greater bandwidth than the sweep. This approach also involvesretention of noise components that reduce the quality of thesubsequently processed data. Finally, Barr discloses phase encoding thesweeps for multiple vibrators, and the use of a different set ofseparation filters for each sweep before stacking the outputs.

U.K. Patent 2,359,363 to Jeffryes (2000) restates the Sallasdisclosures, but with the addition of a filter to remove harmonics fromthe data and from measured vibrator signatures. As noted, filters areundesirable, as they inadequately remove harmonics and othernon-linearities and thus reduce the quality of subsequently processeddata.

There is a need for a method whereby seismic vibrator data can beacquired and processed in a manner that accurately represents the data,which would be derived from an impulse source. The method should involveuse of a deterministic deconvolution that derives from measured vibratormotions. The method should not require addition of white noise tostabilize the processing computations. The method should not require useof a post-processing bandpass filter to eliminate harmonics and noise.The method should retain the correct phase of the underlying data toensure subsequent processing techniques produce accurate results. Themethod should be applicable to arrays of more than one vibrator andprovide for the separation of data recorded from multiple vibratorssimultaneously. The present invention addresses these needs.

SUMMARY

The present invention is a method of processing seismic data in whichone or more seismic vibrators are activated with one or more pilotsignals and vibrator motions are recorded along with seismic data.Vibrator signatures are computed from measured vibrator motions, such asthe ground force signal. A desired impulse response is specified fromeither a measured vibrator motion or from test data or field data from alocation near the location from which the seismic data was acquired. Adeconvolution filter is computed from the impulse response and thevibrator signature. Alternatively, a single separation and deconvolutionfilter is derived from the impulse response and from vibrator signaturesfrom multiple vibrators and sweeps. The deconvolution or deconvolutionand separation filter is used to process the seismic data. The vibratorsare then moved to a new location, and the activation is repeated. Theimpulse response is determined based on an iterative process in whichtime and frequency domain characteristics of it and a sample filter areanalyzed. After application of a separation and deconvolution filter,data are obtained that are suitable for correcting intra-array effectsand for improved noise suppression.

BRIEF DESCRIPTION OF THE DRAWINGS

The features of the present invention will become more apparent from thefollowing description in which reference is made to the drawingsappended hereto.

FIG. 1 depicts a typical data acquisition method that may be used toacquire data for use in embodiments of the present invention.

FIG. 2 depicts a flow chart of the steps for implementation of anembodiment of the present invention.

FIG. 3 depicts a typical pilot sweep signal, ground force signal, andvibrator signature signal where FIG. 3A depicts the entire eight secondrecord of each signal and FIG. 3B depicts an expanded view of the firstone second of each signal.

FIG. 4 depicts a flow chart of the steps for implementation of a secondembodiment of the present invention.

FIG. 5 depicts signal and filter characteristics related to the priorart method of Trantham.

FIG. 6 depicts signal and filter characteristics related to the priorart method of Anderson.

FIG. 7 depicts signal and filter characteristics related to anembodiment of the method of the present invention.

FIG. 8 depicts the inversion operator used in the prior art method ofAllen.

FIG. 9 depicts processed seismic data results for the prior art methodof FIG. 8.

FIG. 10 depicts seismic data results deriving from an embodiment of themethod of the present invention.

FIG. 11 depicts first arrival data recorded for geophones in a well andvarious seismic preprocessing methods, wherein FIG. 11 a showsconventionally correlated data, FIG. 11 b shows data according to themethod of Allen, FIG. 11 c shows the results from an embodiment of thepresent invention, and FIG. 11 d shows impulsive-source data.

FIG. 12 depicts CMP gather data where FIG. 12 a depicts a raw CDP gatherusing four vibrators per station. The four neighboring traces would beautomatically summed by source arrays according to a prior art method.FIG. 12 b depicts the results of applying an embodiment of the presentinvention to individual vibrators and applying static corrections andNMO to each vibrator record.

FIG. 13 depicts the reflector alignments for the two examples of FIG. 12after summing the four neighboring traces to form a source array whereinFIG. 13 a depicts reflector mis-alignment for the prior art method ofFIG. 12 a and FIG. 13 b depicts reflector alignment for the embodimentof FIG. 12 b.

FIG. 14 depicts a comparison of noise suppression results wherein FIG.14 a depicts the results of noise suppression wherein data from fourvibrators at a station are processed according to a prior art method,and FIG. 14 b shows data processed in a supergather according to anembodiment of the present invention in which individual source-receiveroffsets are used.

FIG. 15 depicts the results of processing the data for the two examplesof FIG. 14 in which FIG. 15 a shows obscured reflectors of the prior artmethod of FIG. 14 a and FIG. 15 b shows the visible reflectors accordingto the embodiment of the present invention of FIG. 14 b.

Changes and modifications in the specifically described embodiments canbe carried out without departing from the scope of the invention, whichis intended to be limited only by the scope of the appended claims.

DESCRIPTION OF THE INVENTION

The current invention is a method for improving the quality of vibratorysource data. More specifically, the present invention is a method toacquire and process vibrator seismic data which involves performing avibrator signature deconvolution using vibrator signatures derived frommeasured vibrator motions. Deterministic deconvolution is used insteadof the traditional process of cross correlation and a measured vibratorsignature is used instead of the pilot signal.

The deconvolution filters derived according to the method of the presentinvention remove the vibrator signature from the data, includingharmonics and nonlinear noise, and replace that signature with thedesired impulse response. The filters are designed so as not to amplifynoise for frequencies not included in the vibrator sweep. In addition,the phase distortion caused by the conventional process of adding noiseto stabilize the inversion calculation is reduced.

The method ensures that the output of the initial processing of thevibrator data is similar to that which would be generated by an impulsesource. That initial processing compresses the long wavetrainsassociated with the vibrator sweep into a short-duration impulseresponse with desired phase and amplitude characteristics, therebyallowing first break times to be picked reliably. The present methodensures that the phase of the vibrator data is optimized for subsequentprocessing steps, thereby eliminating the need for subsequent phasecorrections and filtering, and improving the signal-to-noise quality ofthe data. The method uses source signatures for each vibrator and eachsweep to improve the fidelity and resolution of the data, and to allowthe response of each vibrator to be separated. Separating the responseof each vibrator allows the vibrators to be treated as unique sourcepoints, which facilitates intra-array static corrections and noisesuppression.

The method of the present invention builds on prior work in the field ofvibrator data processing, but incorporates characteristics notpreviously recognized to improve the quality of the deconvolved seismicdata. As noted above, for example, Trantham disclosed the deconvolutionto a minimum phase wavelet with a spectrum that goes to zero at bothhigh and low frequencies to match traditional sweeps, and Andersondisclosed the characteristic that the shape of the sweep can be used toreduce the length of the pulse.

A characteristic of the present method is that the impulse responseamplitude spectrum tends to zero faster than does the amplitude spectrumof the vibrator signature—at both high and low frequencies—to preventthe deconvolution filter from becoming unstable in the spectral divisionstep. Prior art methods avoided that instability by adding noise to thesignal and by using an appropriately designed bandpass filter. However,the addition of noise distorts the phase and adds a precursor to thedeconvolved data, and bandpass filters may be inadequate to eliminatenoise, and must be carefully designed to prevent phase distortion. Themethod of the present invention avoids that distortion, since theimpulse response is designed to be rigorously minimum phase withoutneeding a correction, thus ensuring that subsequent processing hascorrect phase. In addition, prior methods, such as that of Barr,intentionally used a desired wavelet with a larger bandwidth than thesweep, which also works against the goal of retaining correct phase.Benefits of the present method include a better ability to pick firstarrival times, and thereby improve the results from checkshot surveys,static correction methods, and tomographic methods. Seismic imagesobtained from the present method have optimum phase and facilitateexcellent ties with well data.

In a multiple vibrator embodiment of the present invention, a fullsolution involving all components of the matrices can be generated usingone filter for all sweeps in a fully coupled derivation. Prior artmethods, such as Barr, use a filter for one sweep, then separate thedata for that sweep. In these prior art methods, separate filters aredesigned to separate subsequent sweeps and separated data from the samelocation are stacked. This approach is equivalent to the use of only thediagonal of the matrix solution to the multivibrator application,results in a loss of data quality, and does not allow the elimination ofnoise or unwanted signal components. Embodiments of the presentinvention avoid these limitations via the full matrix solution.

In a third embodiment, the method of the present invention facilitatesprocessing to be applied for separating the data from each vibrator intoindividual records. Static corrections and differential normal moveout(NMO) can then be applied to each source location. In addition,supergathers can be constructed which improve the capabilities of noiseseparation techniques. The data from each source location can be summed.Alternatively, data can be binned at smaller common depth pointintervals prior to migration, thereby further improve imaging andfocusing.

Embodiments of the present invention will be discussed in the followingin association with the system diagram of FIG. 1, which depicts thegeometry of the data gathering system associated with the presentinvention. In FIG. 1, vibrators 18, 20, 22, and 24 are located, in anon-land application, on trucks 34, 36, 38, and 40, respectively. Signalsthat are generated in the earth by vibrators 18, 20, 22, and 24 arereflected off interface 26 between subsurface layers IM1 and IM2 atpoints 1, 2, 3 and 4, respectively, and are received by receivers 42,44, 46, and 48.

A first embodiment of the method of the present invention, as depictedin the flow chart of FIG. 2, begins with the activation of one or morevibrators with a sweep, step 1. Next, vibrator motions, step 2 andseismic data, step 3 are both recorded. Third, a vibrator signature iscomputed from the measured vibrator motions, step 4. Next, a desiredimpulse response is computed. This involves the specification of adesired impulse response amplitude spectrum, step 5 which is preferablyless than or equal to the vibrator signature spectrum at allfrequencies. Preferably, the amplitudes of the impulse response and thevibrator signature that fall below a small threshold value are set equalto the threshold value, step 6. Alternatively, a small amount ofprewhitening noise can be added at all frequencies (Not depicted in FIG.1). Next, a desired phase for the impulse response is calculated, step7. A deconvolution filter is computed by spectrally dividing thecalculated impulse response by the vibrator signature, step 8. Aniteration 11 is performed if necessary to optimize the characteristicsof the impulse response and the deconvolution filter, step 9. Ifiteration is necessary, step 5 is performed again as well as the stepsafter step 5. The deconvolution filter can then be used for the process13 to deconvolve the received seismic signals, step 10. A characteristicof the present method is that the steps related to the computation ofthe impulse response and the deconvolution filter, once optimized, onlyhave to be performed once and do not have to be repeated for thesubsequent operation of the vibrators at new locations.

As will be understood to those skilled in the art, the impulse responseis specified based on the objectives of the analysis to be performed.For example, for a seismic analysis focusing on a shallow target, itwill be preferable to use a sharper impulse response with higherfrequencies. On the other hand, a deeper target of interest maypreferably require a thicker impulse response with lower frequencies,which penetrate deeper into the earth.

A first embodiment of the present invention will now be described inmore detail. Initially, one or more vibratory sources are used to recordeither a land or a marine seismic survey, and the signals are recordedby one or more detectors, as indicated by the geometry of FIG. 1. Thedetectors can be at the surface as is the case for typical land surveys,suspended within water or at the water bottom as is the case for marinesurveys, or down a wellbore in the case for VSP surveys. Each vibratoris preferably excited, FIG. 2, step 1, with a sweep in which thefrequency of the signal increases or decreases linearly in time;however, any linear, non-linear, or random sweep may be used.Preferably, the sweep should have a low frequency limit at 8 Hz orlower. As will be understood to those skilled in the art, the lowerfrequencies allow for the design of minimum phase impulse responses withrelatively large first lobes. The high frequency limit should be chosenafter tests of different sweep bandwidths and an evaluation of targetreflectors, but will generally be 60 Hz or higher. The sweep signal thatis used to excite the vibrator is also referred to as the pilot signal.

Vibrator motion measurements, FIG. 2, step 2, are made for each sweepand for each vibrator. A preferred method is to record the vibratorground force signal, for example using a mass-weighted sum ofaccelerometers on the vibrator's reaction mass and baseplate or on thevibrator's baseplate stilt assembly. It will be understood that theground force signal is the signal that is used in a feedback mode tocontrol the vibrator excitation. It will also be understood that otheraccelerometer signals or force measurements could be used within thescope of the present invention. FIG. 3A shows an example 8-second pilotsignal, a measured ground force signal, and a vibrator signature signal.The measured ground force is different from the pilot signal because ofharmonics and coupling effects, as clearly indicated in the expandedview, FIG. 3B, of the first 1.0 seconds 31 of the respective signals.

Next, FIG. 2, step 4, a vibrator signature is derived from each measuredmotion. In a preferred method, the signature is computed from the timederivative of the ground force signal in the frequency domain. It isunderstood in the art that the signal imparted into the ground isrelated to the ground force signal. More specifically, the seismic datarecorded by a conventional geophone is in phase with the time derivativeof the ground force. Thus, the vibrator signature can be computed fromthe relationship S(ω)=iωG(ω), where S(ω) is the Fourier transform of thevibrator signature, G(ω) is the Fourier transform of the measured groundforce signal, and the multiplier iω indicates the result in thefrequency domain of the time derivative operation. If more than onevibrator is used, the vibrator signature can be derived from the averageof all of the measured vibrator motions.

It is further understood that the pilot signal is an approximate to theground force signal because of the feedback control. When measurementsare vibrator motions are missing, the reference can be used in place ofthe ground force signal for the design of the vibrator signature. Thiswill allow the same deconvolution to be performed, although harmonicswill not be removed by the process and must be handled separately. Thereference will have zero values at the low and high frequencies.Therefore, special care needs to be taken with the design of the impulseresponse, according to the present invention.

It is a realization of the present invention that it is important forthe vibrator signature to match as closely as possible to the signalthat is actually put into the ground. Use of the time derivative of theground force signal results in less noise and artifacts and less phasedistortion then using the ground force signal itself. Allen and Sallasused the ground force signal for inversion and recovered the timederivative through the process of statistical minimum phasedeconvolution. They, however, required the use of a model trace forphase correction after final processing. With the present invention nosuch correction is required. It is within the scope of the presentinvention to use other modifications, such as the application of asimulated earth filter or attenuation or use of other forcemeasurements, to improve estimates of the vibrator signature.

In the vibrator signature computation the product iωG(ω) may also bemultiplied by a scale factor C to normalize the amplitude of thevibrator signature amplitude spectrum to unity. In addition, vibratorsignature values less than a threshold value T may be set equal to thatthreshold T, such that S(ω) may, for example, be expressed as follows:S(ω)=iωCG(ω) when G(ω)>T/C  (1)S(ω)=iT/C when G(ω)≦T/C  (2)Threshold T may be specified in absolute value terms, or as a percent ofthe peak value, or otherwise.

The desired impulse response is constructed next. In the frequencydomain this impulse response will have the formI(ω)=A ₁(ω)e ^(−i(φ) ¹ ⁾  (3)where A₁(ω) and φ₁(ω) are the amplitude and phase of the desiredimpulse, respectively. Either a measured vibrator motion from the surveybeing performed may be used to construct the impulse response, or theimpulse response may be derived before the survey based on test data oron data from a similar location. First, FIG. 2, step 5, an impulseamplitude spectrum A₁(ω) that is appropriate for the bandwidth of thesweep signal is constructed. It will be understood that the amplitudespectrum should preferably have a peak in the middle of the sweep bandand should not have sharp edges. The amplitude spectrum of Anderson ispreferable. However, the spectrum of Trantham may also be used, as mayother spectral forms that will be known to those skilled in the art. Forconvenience in the following, the amplitude spectrum will be assumed tohave a peak amplitude of one, with the amplitude decreasing at higherand lower frequencies relative to the center portion of the bandwidth.

One implication of designing the amplitude spectrum based on thevibrator signature and not on the sweep signal is that the peak of theimpulse amplitude spectrum may be at a different frequency from the peakof the sweep amplitude. For example, when a time derivative is employedto compute the signature, as described above, higher frequencies areboosted, thus allowing the use of a higher frequency impulse. In suchcases, however, the amplitudes at the lower frequencies often will bereduced compared to that in a standard linear sweep, in order to make ashorter duration pulse. In addition, it is well-known, for example asdiscussed by Anderson and Trantham, that a square impulse responsecorresponds to a longer pulse, and tapering or shaping shortens thepulse. It is also understood that harmonics in the ground force cancontribute to the tails of the frequency response, but do not serve tosubstantially increase the frequency or bandwidth in the desired impulseresponse. It is preferable that the amplitude spectrum for the desiredimpulse at the low and high frequencies is modified by smoothly taperingthe amplitude to the threshold value at both low and high frequencies.This ensures that the filter amplitudes remain near or below unity (zeroamplification) outside the sweep band of the vibrator.

Next, FIG. 2, step 6, a threshold value is selected for the impulseresponse amplitude spectrum, such that all values of the amplitudespectrum that fall below the threshold value are set equal to thatvalue. Alternatively, prewhitening can be performed in which a smallconstant amplitude, essentially white noise, is added for allfrequencies.

One reason for use of a threshold value in step 6 is to ensure accuratecalculation of the impulse response phase spectrum, φ₁(ω), FIG. 2, step7. In this step both a minimum phase and zero phase impulse are computedfrom the amplitude spectrum. The minimum phase impulse is desirable forsuch analyses as first arrival determination. The zero phase impulsewill represent that of the data after final processing and is thereforepreferable for data interpretation. The zero phase impulse is computedby taking the inverse Fourier transform of the amplitude spectrum, andthat computation assumes that the phase is zero throughout thebandwidth. The minimum phase impulse is computed by taking the Hilberttransform of the natural logarithm of the amplitude spectrum. Becausethe natural logarithm calculation requires a non-zero impulse responseamplitude spectra value at all frequencies, the threshold value appliedin step 6 ensures the non-zero amplitude condition throughout thebandwidth of the impulse response spectrum. The inverse Fouriertransform is then used to compute the minimum phase impulse response.

The amplitude spectrum of the deconvolution filter, FIG. 2, step 8, isconstructed by dividing the amplitude spectrum of the impulse responseby the amplitude spectrum of the vibrator signature. A scale factor maybe employed for the amplitude spectrum of the vibrator signature toensure that the amplitude of the impulse spectrum matches that of thevibrator signature in the region of the center frequency of thebandwidth. Alternatively, the root mean square amplitude in a specifiedfrequency window may be used to match the amplitudes. The amplitudes atthe low and high frequency ends of the amplitude spectrum of thevibrator signature that fall below a threshold value are set equal to athreshold to facilitate subsequent use of the deconvolution filter. Thethreshold is generally set at a percent of the peak value of thespectrum, and should be equal to or larger than the threshold used forthe impulse response. This threshold only affects the low and highfrequencies in the subsequent processing, and is used for analyticstability.

In step 9 of FIG. 2, both the impulse response and the deconvolutionfilter are studied to determine whether their characteristics aresuitable for the intended analysis. The time impulse response iscomputed by inverse Fourier transform, and is evaluated as tosuitability for the analysis to be performed. The impulse time curvesshould have few lobes, be short in time, and the minimum phase impulseshould have most of its energy in the early part of the pulse. Theprocess is iterated until optimum impulse response spectra are obtainedin terms of both the time response of the impulses and the amplificationof the filter. It will be understood to those skilled in the art thatthis iteration process will involve evaluation of the time domaincharacteristics of the impulse pulses, both zero phase and minimumphase, and evaluation of the frequency domain characteristics of theimpulse response amplitude spectrum, in particular in the high and lowfrequency portions of the spectrum.

The iteration will also include an evaluation of the frequency domaincharacteristics of the amplitude spectrum of the filter, again withparticular emphasis on the high and low frequency portions of thespectrum. Finally, the iteration may involve selective application ofthe filter to actual seismic data to determine the filter's ability toeliminate noise in the processed data and the tendency of the processeddata to show the presence of precursors. It will be understood thatamong the characteristics that will be studied in an actual seismic datatest of the filter will be the ability to determine first arrival times,and the extent to which the shape of the wavelet is reasonably clear andclean. As noted above, among the general traits that the impulseamplitude spectrum will typically contain is an amplitude which is lessthan that of, and which trends to a zero magnitude faster than does, thevibrator signature spectrum in the low frequency range.

The impulse response will also typically trend to zero magnitude athigher frequencies at a rate sufficient that the response's magnitude isless than, and preferably set to the desired threshold value, themagnitude of the vibrator signature in the frequency range whereharmonics begin to be observed. The filter resulting from the impulseresponse will preferably have amplitudes less than one in the lowfrequency range, such as below 8 Hz, and in the high frequency end abovethe highest frequencies in the sweep. Filters. according to the presentinvention will also have magnitudes in the high frequency range smallerthan the amplitudes of filters using prior art methods.

In the frequency domain, the deconvolution filter F(ω) is the desiredimpulse response in the frequency domain I(ω) divided by the vibratorsignature S(ω), or $\begin{matrix}{{F(\omega)} = {\left( \frac{A_{I}(\omega)}{\omega\quad{A_{g}(\omega)}} \right)\left( {\mathbb{e}}^{- {i{({\phi_{I} - {\pi/2} - \varphi_{g}})}}} \right)}} & (4)\end{matrix}$where A₁(ω) and φ₁(ω) are the amplitude and phase of the desired impulseand A_(G)(ω) and φ_(G)(ω) are the amplitude and phase of the vibratorground force signal. The time derivative computation for the vibratorsignature is the amplitude scaling by 1/ω and the phase shift by 90degrees or π/2.

The deconvolution filter is then applied to process the seismic datafrom step 3 in the frequency domain, FIG. 2, step 10. It is importantfor minimum phase causal results that for every amplitude component ofthe deconvolution filter the equivalent minimum phase function is alsoapplied. Conversely for every phase correction, such as a 90-degreephase rotation, the corresponding amplitude correction, such as 1/ω,must also be applied. As described above, the present inventionsatisfies both of these criteria.

In a second embodiment of the present invention, a vibrator signaturedeconvolution is performed within a matrix separation scheme for anumber of vibrators operated simultaneously. In this embodiment, aplurality of vibrators is then employed to perform a number of sweeps.The number of sweeps should be equal to or greater than the number ofvibrators. Measurements are made of both the motions of the vibratorsand the received seismic signals. Next, vibrator signatures arecomputed, with the amplitudes clipped at a minimum threshold value. Adesired impulse response is constructed next, as described above. Adeconvolution matrix operator is generated in the frequency domain thatseparates the earth response for each vibrator and replaces theindividual vibrator signatures with the desired impulse response. Thatoperator is applied to deconvolve the seismic data and to separate thatdata according to the individual vibrator locations.

In this second embodiment, multiple sweeps, in a number equal to orgreater than the number of vibrators, are simultaneously obtained fromall vibrators, and a matrix separation and vibrator deconvolution methodis employed. Preferably, the separation method involves use of phaseencoding. For example, with three vibrators and three sweeps, a sequencemay be constructed where on each sweep one vibrator is operated 90degrees out of phase with the other vibrators, such as in the following:Sweep Vibrator 1 Vibrator 2 Vibrator 3 1 90 0 0 2 0 90 0 3 0 0 90

It would also be possible to have a fourth sweep in which all vibratorstime and at the same phase during one sweep, such as: Sweep Vibrator 1Vibrator 2 Vibrator 3 1 90 0 0 2 0 90 0 3 0 0 90 4 0 0 0

The above phase encoding may also be superimposed along with variphasingto further reduce harmonics in the case where the vibrator signaturedoes not perfectly match the harmonics put into the ground. Variphasinginvolves phase rotations of a factor of 2π/M where M is the number ofsweeps, as described by E. Rietsch in “Reduction of Harmonic Distortionin Vibratory Source Records,” Geophysical Prospecting, v. 29, pp.178-188, 1981. This suppresses all harmonics up to and including orderM. For use in this application the variphase angle is summed with the 90degrees phase shift. For example for M=4 sweeps, the variphase rotationsare 0, 90, 180, and 270 degrees. It is understood that the phaserotations can be performed in any order, and that higher level multiplesof 2π in the above factor may be employed within the scope of thepresent invention. For example, adding a ninety degrees phase encodingfor one vibrator at a time yields the following. The first sweep is aphase angle of 0 degrees for all vibrators except the first vibratorwhich sweeps at 0+90 or 90 degrees. The second sweep is at a phase of180 for all vibrators except for the second vibrator that sweeps at180+90 or 270 degrees, and so forth for the other two sweeps. SweepVibrator 1 Vibrator 2 Vibrator 3 Vibrator 4 1 90 0 0 0 2 90 180 90 90 3180 180 270 180 4 270 270 270 0It will be understood that other encoding methods may also be employed.

This second embodiment will now be discussed in more detail, inassociation with FIG. 4. Consider N vibrators radiating M≧N sweeps intothe earth, FIG. 4, step 50) and resulting in M recorded data traces(FIG. 4, step 54). It is desired to obtain, by solving a set of linearequations that finds the set of N earth reflectivities, the operatorthat best predicts the recorded data, based on the known MN vibratorsignatures. Referring to FIG. 1, which shows four vibrators 18, 20, 22,24 radiating a different signature s into the ground, s₁, s₂, s₃, s₄.Each signature is convolved with a different earth reflectivity sequencee₁, e₂, e₃, e₄ (For example, as resulting from the reflections at points1, 2, 3, and 4 in FIG. 1). The earth sequence can include reflectors,multiples, and near-surface effects.

Traces d_(i)(t) recorded at each geophone 42, 44, 46, 48 are a sum ofthe signature-filtered earth reflectivities under each vibrator. Thedata trace d_(i)(t) recorded for sweep i is: $\begin{matrix}{{d_{i}(t)} = {\sum\limits_{j = 1}^{N}{{s_{ij}(t)} \otimes {e_{j}(t)}}}} & (5)\end{matrix}$where s_(ij)(t)=sweep i from vibrator j and e_(j)(t)=the earthreflectivity seen by vibrator j.

In the frequency domain, this expression becomes: $\begin{matrix}{{D_{i}(\omega)} = {\sum\limits_{j = 1}^{N}{{S_{ij}(\omega)}{E_{j}(\omega)}}}} & (6)\end{matrix}$which, in matrix form for M sweeps and N vibrators, can be expressed$\begin{matrix}{{{\begin{bmatrix}S_{11} & S_{12} & \vdots & S_{1N} \\S_{21} & S_{22} & \vdots & S_{2N} \\S_{31} & S_{32} & \vdots & S_{3N} \\S_{41} & S_{42} & \vdots & S_{4N} \\\vdots & \vdots & \vdots & \vdots \\S_{M\quad 1} & S_{M\quad 2} & \vdots & S_{MN}\end{bmatrix}\begin{bmatrix}E_{1} \\E_{2} \\\vdots \\E_{N}\end{bmatrix}} = \begin{bmatrix}D_{1} \\D_{2} \\D_{3} \\D_{4} \\\vdots \\D_{M}\end{bmatrix}}{or}} & (7) \\{{S\quad\overset{\rho}{E}} = \overset{\rho}{D}} & (8)\end{matrix}$

If the number of sweeps is equal to the number of vibrators, the earthresponse can be specified in terms of a filter F: $\begin{matrix}{\overset{\rho}{E} = {{F\overset{\rho}{D}\quad{where}\quad F} = S^{- 1}}} & (9)\end{matrix}$Filter F is the inverse of the matrix of vibrator signatures. In thisprior art solution, inversion filter F separates the response of eachvibrator and compresses the sweeps to impulses. However, as discussedabove in association with a first embodiment of the present invention,it more desirable to use a deconvolution filter which incorporates aminimum phase impulse response I to determine the earth response. Inmatrix form, a preferred filter according to the present invention isF=I(S ⁻¹)  (10)

Note that equation (1) is the matrix equivalent of equation (4)discussed above. With reference to FIG. 4, the vibrator motions arerecorded, step 52, and vibrator signatures computed, step 56, in amanner analogous to that described above in association with FIG. 2. InFIG. 4, step 58, the impulse response is constructed and the matrix ofvibrator signatures are inverted 60 according to the procedure describedabove in association with FIG. 2, steps 5, 6, and 7. That impulseresponse is used in equation (10) to determine the deconvolution filterF, step 62, and to thereafter separate and deconvolve the recordedseismic data, step 64. It will be understood to those skilled in the artthat the iteration step of FIG. 1, step 9, which is not depicted in FIG.4, may be included in multi-vibrator embodiments of the presentinvention, as desired or necessary. If there are more sweeps thanvibrators, the multi-vibrator problem is overdetermined, and a solutionmust be determined by a least squares analysis. The normal equations are$\begin{matrix}{{{S^{*}S\quad\overset{\rho}{E}} = {S^{*}\overset{\rho}{D}}}{\overset{\rho}{E} = {\left( {S^{*}S} \right)^{- 1}S^{*}\overset{\rho}{D}}}} & (11)\end{matrix}$where S* is the conjugate transpose of the vibrator signature matrix S.Equations 11 involve a prior art method that follows the disclosure ofSallas et al. in U.S. Pat. No. 5,721,710. According to the presentinvention however, it is preferable to use a deconvolution to remove thevibrator signature S_(ij) and replace it with an impulse response I. Inthis embodiment the deconvolution filter F becomesF=(S*S)⁻¹(S*I)  (12)

It will be understood to those skilled in the art that there are anumber of ways to solve for filter F in equation (12). One such methoduses a LU-decomposition of the matrix (S*S) and forward and backsubstitution to find the operator (S*S)⁻¹(S*I). Further details on thisand other such methods may be found in Numerical Recipes, W.H. Press etal, Cambridge University Press, 1986. Once the filter matrix F isgenerated, it is applied to each data trace in the frequency domain toobtain a separated record for each vibratory source.

After the data are separated and inverted according to the method of thepresent invention, various processes can be applied. The method of thepresent invention uses the best estimate of the vibrator signature todeconvolve the data, but it is recognized that the actual energy putinto the ground may be related to the computed signature by otherminimum phase processes. For example, the actual signature may beaffected by near-surface effects. The process of surface-consistentdeconvolution or spiking deconvolution can remove these effects, asdiscussed by Allen. Subsequent processes which may also be appliedinclude, but are not limited to, noise suppression, divergencecorrection, and static corrections.

In multi-vibrator applications of the present invention, separation ofthe data into unique source points for each vibrator permits specialprocessing techniques to improve data quality. Conventionally, forexample, a number of vibrators are operated simultaneously at a sourcestation to form a source array. The array may suppress some ground-rollnoise, but because the number of vibrators is usually small, noisesuppression by the array is typically not very effective. Also, if thevibrators are located at different elevations, the reflections mayarrive at receivers with small time differences.

The result of the time differences will limit the higher frequencycomponents of the processed data. Finally, the trace spacing of theprocessed data is limited by the source and receiver station interval.By separating the data into unique source points according toembodiments of the present invention, the reflectors can be alignedbefore summing the data into arrays. Source static corrections for eachvibrator can then be computed by any conventional method and applied tothe data. In addition, correction can be made for differential normalmoveout between the vibrator units because of slightly different offsetsfor each source. In addition, source generated noise can then besuppressed by using the unique offsets for each source.

The data recorded by a set of vibrators may be sorted into a supergatherby sorting by the unique offset for each source and each receiver. Theextra spatial sampling in the supergather will permit coherent noise tobe removed by processes such as FK filtering that were aliased at theoriginal spacing. After removing noise and aligning the reflectors,array forming can be performed. Alternatively, the data can be binnedinto finer CDP bins using the separated source interval and input intoprestack migration. For example, if the original source and receiverinterval is 50 m, then the CDP interval is 25 m. However, if vibratorsare used per source interval, and the records separated, then the datacan be binned at 6.25-m intervals. The output can be either the finertrace spacing or a coarse spacing, whichever is desired.

The benefits of designing deconvolution filters according to embodimentsof the present invention can be demonstrated using the signals of FIG.3. Recall, that FIG. 3 shows a pilot sweep signal and a ground forcesignal from a standard 8-second linear sweep. The frequency range isfrom 8 to 128 Hz. As can be seen, the ground force is not as smooth andcontinuous as is the pilot because of harmonics. For reference, FIG. 3also shows a vibrator signature, computed using the derivative of theground force.

FIG. 5 illustrates the design of a filter 506 from a vibrator signature504 using the impulse amplitude spectrum 502 from Trantham. An advantageto this amplitude spectrum is that its magnitude goes to zero at a highand low frequency, as does the sweep. However, note that the impulseamplitude spectrum 502 is not a good match to that of the vibratorsignature 504, which causes a peak in the amplitude of filter 506 at arelatively low frequency compared to the center of the bandwidth. Inaddition, the impulse responses 510 and 512 corresponding to the impulseamplitude spectrum 502 exhibit some ringing at longer times, and it willbe understood in the art that such ringing characteristics are notdesired.

In FIG. 6, an amplitude spectrum 602 from Anderson is used. Thisspectrum does not go to zero when the sweep goes to zero, and note thatthe amplitude spectra better matches that of the signature 704. However,it will be observed that filter 606 has substantially high amplitudepeaks at both the low and high frequency ranges, with filter amplitudesmuch greater than 1.0. As a result the filter amplifies background noiserecorded by geophones, 608. Finally, the impulse responses 610 and 612are improved as compared to Trantham.

FIG. 7 depicts a filter 706 and impulse responses 710 and 712 derivingfrom an embodiment of the present invention. In this solution a taper toa minimum threshold value was added at both low and high frequency toamplitude spectrum 702 (The minimum threshold value is not visible inFIG. 7). In addition, the taper at the upper portion of the frequencyband was altered such that spectrum 702 is below the amplitude ofsignature 704 for essentially the entire upper portion of the frequencyrange. The result is a filter 706 with substantially reduced amplitudesat both the upper and lower portions of the frequency range, and cleanthree-lobed impulse wavelets 710 and 712 each having a large primarylobe.

FIGS. 8-15 illustrate the advantages of use of embodiments of thepresent invention. FIG. 8 shows a prior art inversion designed from thesweep in FIG. 3 according to the patent to Allen. Because the sweepstarts at 8 Hz and stops at 128 Hz 81, the inverse of the sweep 83results in large gains at frequencies below 8 Hz and above 128 Hz. Theinverse of the sweep and noise 85 is also shown. The result of applyingthis inversion operator to seismic data is shown in FIG. 9. FIG. 9 a,which depicts data recorded by geophones in a well, exhibit firstarrivals that are obscured by the low frequency noise that is amplifiedby the inversion filter. This is also evident in the spectrum in FIG. 9b, which also shows large magnitudes at low frequency. In comparison,FIG. 10 shows the result of application of a deconvolution filteraccording to an embodiment of the present method to the same geophonewell data. The data in FIG. 10 a are clean and noise-free, clearlyshowing first-arrivals and reflections. The spectrum in FIG. 10 b showsminimal signal energy outside the sweep band of 8-128 Hz. Thus, there isno need for subsequent filtering to remove noise at the low and highfrequencies, an improvement of the present invention over the prior art.

As discussed above, the conventional method for stabilizing an inversionfilter is to prewhiten the signature by adding white noise, resulting inan inversion filter as shown in FIG. 8. However, with prewhitening, anamplitude correction is made to the filter without a corresponding phasecorrection. Although for small added noise amplitudes, less than 1percent, the effect may be small, higher noise levels, in the range of 3to 5 percent, may be needed to prevent large gains at the lower andhigher frequencies. The result is phase distortion and precursors in theprocessed data. It will be understood in the art that there should be noprecursor for first arrivals with truly minimum phase data. FIG. 11depicts first arrival data for various preprocessing method, and furtherdemonstrates the benefits of the present invention. FIG. 11 d shows theresults from impulsive-source data recorded by geophones in a well. FIG.11 a shows conventionally correlated data, and the poor match to theimpulsive-source data is clearly evident. FIG. 11 b involved aninversion according to the method of Allen, and although the match tothe impulsive-source data is improved the result is still noisy andshows a precursor. FIG. 11 c shows the results of applying adeconvolution filter according to an embodiment of the presentinvention. Precursor energy is small and the wavelets are a good matchto the impulsive-source data.

FIGS. 12-15 illustrate the improved seismic data processing resultswhich derive from use of embodiments of the method of the presentinvention. FIGS. 12 and 13 demonstrate the advantages of the ability ofembodiments of the present invention for use of unique source pointsafter separation of the individual recordings from each vibrator. FIG.12 a shows a raw CDP gather using 4 vibrators per station. Withconventional seismic data recording, the vibrators form a source arrayand data from each group of 4 traces are combined in the field as theyare recorded. In processing, conventional static corrections and NMOmoveout are applied. The resulting CMP gather in FIG. 13 a clearly showsthat the reflectors are not aligned. However, by separating theindividual vibrators using an embodiment of the present invention, withstatic corrections and NMO applied to each separated vibrator record, asshown in FIG. 12 b, the data which results, after forming the sourcearray by summing the 4 traces depicted in the gather in FIG. 13 b, showsmuch better reflector continuity as compared to FIG. 13 a.

In addition, the present invention facilitates better noise suppression,as is illustrated in FIGS. 14 and 15. FIG. 14 a shows the results of aprior art method involving 4 vibrators per station, and a large amountof aliased ground-roll noise 141 is present. When the data is summed,noise suppression techniques cannot remove much of the aliased groundroll without removing signal. The result is shown in FIG. 15 a, in whichreflectors are obscured by noise. In contrast, FIG. 14 b shows arrangingthe data in a supergather using individual source-receiver offsetsaccording to an embodiment of the method of the present invention. Muchof the ground roll noise is no longer aliased, and can be removed priorto forming the source array. The result, shown in FIG. 15 b, is data inwhich reflectors can be seen at near offsets, no longer obscured by thenoise.

It should be understood that the preceding is merely a detaileddescription of specific embodiments of this invention. Other embodimentsmay be employed and numerous changes to the disclosed embodiments may bemade in accordance with the disclosure herein without departing from thespirit or scope of the present invention. Furthermore, each of the aboveembodiments is within the scope of the present invention. The precedingdescription, therefore, is not meant to limit the scope of theinvention. Rather, the scope of the invention is to be determined onlyby the appended claims and their equivalents.

1. A method of processing seismic data generated by a seismic vibratorcomprising the steps of (a) measuring a signal related to the vibratormotion; (b) computing a vibrator signature for said measured vibratormotion signal; (c) specifying a desired seismic data processing impulseresponse, wherein the high and low frequency portions of an amplitudespectrum of said impulse response taper to zero at a rate faster thandoes the high and low frequency portions of an amplitude spectrum ofsaid vibrator signature; (d) computing a deconvolution filter from theratio of the desired impulse response and the computed vibratorsignature; (e) processing said seismic data using said deconvolutionfilter.
 2. The method of claim 1, wherein said vibrator motion signal isthe vibrator ground force signal.
 3. The method of claim 1, wherein saidvibrator motion signal is approximated by a vibrator pilot signal. 4.The method of claim 1, wherein small values of an amplitude spectrum ofsaid impulse response are set equal to a threshold value.
 5. The methodof claim 1, wherein small values of an amplitude spectrum of saidimpulse response are set less than a threshold value.
 6. The method ofclaim 1, wherein the vibrator signature is computed from the timederivative of said ground force signal.
 7. The method of claim 1,wherein said computation of said deconvolution filter involves aniteration in which zero phase and minimum phase impulse responses areanalyzed to determine the suitability of said deconvolution filter toprocess high and low frequency data content in said seismic data.
 8. Themethod of claim 1, wherein a plurality of seismic vibrators is used, anda resulting matrix of vibrator signatures is inverted and used todetermine said deconvolution filter.
 9. The method of claim 1, wherein asmall amount of pre-whitening noise is added to all frequency values ofan amplitude spectrum of said impulse response
 10. The method of claim1, wherein a plurality of vibratory sources are used to generate seismicdata and a plurality of detectors are used to record the seismic data,wherein each vibrator is excited with a frequency sweep.
 11. The methodof claim 10 wherein a full solution involving all components of thematrix is generated by using one filter for all sweeps in a fullycoupled derivation of the full matrix solution.
 12. The method of claim1 wherein the frequency sweep increases linearly in time.
 13. The methodof claim 1 wherein the frequency sweep decreases linearly in time. 14.The method of claim 1 wherein the frequency sweep is a nonlinear sweep15. The method of claim 1 wherein the frequency sweep is a random sweep16. The method of claim 10 wherein different vibrators are energized bydifferent sweeps which are phase encoded.
 17. The method of claim 10wherein different vibrators are energized by sweeps in which onevibrator a time is energized by a sweep with a 90 degree phase rotationrelative to the phases of the other vibrators
 18. The method of claim 10wherein multiple sweeps are used and the sweeps include phase rotationsof 360/N degrees where N is an integer.
 19. The method of claim 17wherein multiple sweeps are used and the sweeps include phase rotationsof 360/N degrees where N is an integer.
 20. The method of claim 10wherein the different vibrators are energized by sweeps coveringdifferent frequency ranges at different times.
 21. The method of claim10 wherein the different vibrators are energized by sweeps which startat different times.
 22. The method of claim 1 wherein the location ofthe detectors is selected form the group comprising detectors on thesurface of the earth, detectors suspended in the water, detectors on thewater bottom, detectors in a wellbore, and any combination thereof. 23.The method of claim 10 further comprising: (a) separating the data fromeach vibrator into individual records for each source location; (b)applying static corrections and differential normal moveout (NMO)correction to each source location; (c) summing the data for each sourcelocation.
 24. The method of claim 10 further comprising (a) separatingthe data from each vibrator into individual records for each sourcelocation; (b) constructing supergathers to improve the noise separationtechniques, and (c) summing the data for each source location
 25. Themethod of claim 10 further comprising separating the data into bins atsmall common depth point intervals and migrating the data therebyimproving imaging and focusing.
 26. A method of processing seismic datagenerated by at least two seismic vibrators with a number of sweeps atleast equal to the number of vibrators comprising the steps of (a)measuring a vibrator motion signal for each vibrator; (b) measuring theseismic signal; (c) computing a vibrator signature for said measuredvibrator motion signal; (d) specifying a desired seismic data processingimpulse response, wherein the high and low frequency portions of anamplitude spectrum of said impulse response taper to zero at a ratefaster than does the high and low frequency portions of an amplitudespectrum of said vibrator motion signal; (e) computing a deconvolutionmatrix in the frequency domain that separates the earth response foreach vibrator and replaces the individual vibrator signatures with thedesired impulse response; (f) processing said seismic data andseparating the data according to the individual vibrators using saiddeconvolution matrix.
 27. The method of claim 26, wherein said vibratormotion is the vibrator ground force signal.
 28. The method of claim 26wherein said vibrator motion is approximated by a vibrator pilot signal.29. The method of claim 26, wherein small values of an amplitudespectrum of said impulse response are set equal to a threshold value.30. The method of claim 26, wherein small values of an amplitudespectrum of said impulse response are set less than a threshold value.31. The method of claim 26, wherein the vibrator signature is computedfrom the time derivative of said ground force signal.
 32. The method ofclaim 26, wherein said computation of said deconvolution filter involvesan iteration in which zero phase and minimum phase impulse responses areanalyzed to determine the suitability of said deconvolution filter toprocess high and low frequency data content in said seismic data. 33.The method of claim 26, wherein a small amount of prewhitening noise isadded to all frequency values of an amplitude spectrum of said impulseresponse
 34. The method of claim 26 wherein the frequency sweepincreases linearly in time.
 35. The method of claim 26 wherein thefrequency sweep decreases linearly in time.
 36. The method of claim 26wherein the frequency sweep is a nonlinear sweep.
 37. The method ofclaim 26 wherein the frequency sweep is a random sweep.
 38. The methodof claim 26 wherein different vibrators are energized by sweeps whichare phase encoded.
 39. The method of claim 26 wherein differentvibrators are energized by sweeps in which one vibrator a time isenergized by a sweep with a 90 degree phase rotation relative to thephases of the other vibrators
 40. The method of claim 26 whereinmultiple sweeps are used and the sweeps include in any order of phaserotations of 360/N degrees where N is an integer.
 41. The method ofclaim 39 wherein multiple sweeps are used and the sweeps include in anyorder of phase rotations of 360/N degrees where N is an integer.
 42. Themethod of claim 26 wherein the different vibrators are energized bysweeps covering different frequency ranges at different times.
 43. Themethod of claim 26 wherein the different vibrators are energized bysweeps which start at different times.
 44. The method of claim 26wherein the location of the detectors is selected form the groupcomprising detectors on the surface of the earth, detectors suspended inthe water, detectors on the water bottom, detectors in a wellbore, and acombination thereof.
 45. The method of claim 26 further comprising: (a)separating the data from each vibrator into individual records for eachsource location; (b) applying static corrections and differential normalmove-outs (NMO) to each source location; (c) summing the data for eachsource location.
 45. The method of claim 26 further comprisingconstructing supergathers to improve the noise separation techniques.46. The method of claim 26 further comprising separating the data intobins at small common depth point intervals and migrating the datathereby improving imaging and focusing.